###### CPS paper on trust and risky decision-making ######
###### Replication file #####
## Run in R version 4.0.3 (Bunny-Wunnies Freak Out)
#### Please note that the data and code provided are taken from a much larger project on trust and governance. 
#### If you are interested in finding out more about the rest of the data collected, please contact the author.

#### Directory and Library Prep ####
setwd("INSERT DIRECTORY")
getwd()

install.packages("pacman")
library(pacman)
pacman::p_load(interplot, oddsratio, rmarkdown, psych, stargazer, ggplot2, foreign, ggridges, 
               gridExtra, survival, knitr, broom, purrr, lavaan, Rmisc, car, tidyverse, dplyr, 
               survival, lme4, lmerTest, reshape, DescTools, blorr, margins, lmtest, sandwich, MASS, Hmisc,
               weights, anesrake, survey, data.table, jtools, interactions, QuantPsyc, performance, see, 
               patchwork, ggeffects, ggcorrplot)

## Read data
MPs <- read.csv("CPS_REPLICATION_DATA.csv", stringsAsFactors = FALSE)
names(MPs)

#### Housekeeping ####
## Create a rescaling formula for housekeeping on continuous variables
range01 <- function(x, ...){(x-min(x, ...))/(max(x, ...)-min(x, ...))}

#### RECODE: Ethnicity ####
table(MPs$Ethnicity)
MPs$Ethnicity1 <- MPs$Ethnicity

MPs$Ethnicity1[MPs$Ethnicity1=="Asian South African"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Any other black background"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Any other mixed background"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Other ethnic background"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Black South African"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Black Caribbean"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Chinese"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Indian"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Other ethnic group"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="Maori"] <- "BIME"
MPs$Ethnicity1[MPs$Ethnicity1=="White and Asian"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1=="White Australian"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1=="White British"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1=="Any other white background"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1=="White"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1=="White South African"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1=="Prefer not to say"] <- NA
MPs$Ethnicity1[MPs$Ethnicity1=="White/New Zealand European"] <- "White"
MPs$Ethnicity1[MPs$Ethnicity1==""] <- NA
table(MPs$Ethnicity1)

#### RECODE: Education ####
table(MPs$Education)
MPs$Education1 <- MPs$Education

MPs$Education1[MPs$Education1=="(Modern) Apprenticeship, Advanced (Modern) Apprenticeship, SVQ/NVQ/Key Skills Level 1 and 2 City and Guilds Craft/Inter"] <- "Apprenticeship"
MPs$Education1[MPs$Education1=="3-4 year University/CNAA first Degree (BA, BSc., BEd., BEng. etc)"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="5 year University/CNAA first Degree (MB, BDS, BV etc)"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="City and Guild certificate"] <- "Apprenticeship"
MPs$Education1[MPs$Education1=="Clerical and commercial"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="CSE grade 1, GCE O level, GCSE, School Certificate"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="CSE grades 2-5"] <- "None"
MPs$Education1[MPs$Education1=="Don't know"] <- "None"
MPs$Education1[MPs$Education1=="Edexcel/BTEC/BEC/TEC - Higher National Certificate (HNC) or equivalent"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="Foundation Degree (FdA, FdSc etc)"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="GCE A level or Higher Certificate"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="HE Access"] <- "A-Levels/Vocational Diploma"
MPs$Education1[MPs$Education1=="Masters Degree, M.Phil, Post-Graduate Diplomas and Certificates"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="No formal qualifications"] <- "None"
MPs$Education1[MPs$Education1=="None of these"] <- "None"
MPs$Education1[MPs$Education1=="Nursing certificate, Teacher training, HE Diploma, Edexcel/BTEC/BEC/TEC - Higher National Diploma (HND), OCR/RSA - Highe"] <- "A-Levels/Vocational Diploma"
MPs$Education1[MPs$Education1=="Nursing qualification (eg SEN, SRN, SCM, RGN)"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="ONC"] <- "Apprenticeship"
MPs$Education1[MPs$Education1=="Other"] <- "None"
MPs$Education1[MPs$Education1=="Other technical, professional or higher qualification"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="Ph.D, D.Phil or equivalent"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="Recognised trade apprenticeship completed"] <- "Apprenticeship"
MPs$Education1[MPs$Education1=="Refusal"] <- "None"
MPs$Education1[MPs$Education1=="Scottish Higher Certificate"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="Scottish Ordinary/ Lower Certificate"] <- "None"
MPs$Education1[MPs$Education1=="Teaching qualification (not degree)"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="University diploma"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="University or CNAA first degree (eg BA, B.Sc, B.Ed)"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="University or CNAA higher degree (eg M.Sc, Ph.D)"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="Vocational A-level (AVCE), GCE Applied A-level, NVQ/SVQ Level 3 GNVQ/SNVQ Advanced, Edexcel/BTEC/BEC/TEC (General/Ordina"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="Youth training certificate/skillseekers"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1==""] <- NA
MPs$Education1[MPs$Education1=="City and Guild certificate – advanced"] <- "Apprenticeship"
MPs$Education1[MPs$Education1=="A-Levels/Vocational Diploma"] <- "Secondary (school leavers) Qualification or Vocational Diploma"
MPs$Education1[MPs$Education1=="Bachelor's degree, Advanced Diplomas, Post Graduate Certificate and B-tech"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="College, CEGEP or other non-university certificate or diploma"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="Doctorate"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="Grade 12 (National Senior Certificate) and National (vocational) Cert. level 4"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="High school diploma or equivalency certificate"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="Higher Certificates and Advanced National (vocational) Cert."] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="Honours degree, Post Graduate diploma and Professional Qualifications"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="Level 10 – Doctoral Degree"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="Level 5 – Diploma"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="Level 7 – Bachelor Degree"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="Level 8 – Bachelor Honours Degree, Graduate Certificate, Graduate Diploma"] <- "Higher Education Degree"
MPs$Education1[MPs$Education1=="Level 9 – Masters Degree"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="Master's degree"] <- "Postgraduate Degree"
MPs$Education1[MPs$Education1=="National Diploma and Advanced certificates"] <- "Tertiary Diploma"
MPs$Education1[MPs$Education1=="NZ school certificate in one or more subjects"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="Senior Secondary Certificate of Education"] <- "Secondary (school leavers) Qualification"
MPs$Education1[MPs$Education1=="University certificate or diploma below bachelor level"] <- "Tertiary Diploma"

table(MPs$Education1)

MPs$Education2 <- factor(MPs$Education1, 
                         levels = c("None", 
                                    "Secondary (school leavers) Qualification",
                                    "Tertiary Diploma",
                                    "Higher Education Degree",
                                    "Postgraduate Degree"), ordered = TRUE)

#### RECODE: Ideology ####
describe(MPs$EconIdeology)
describe(MPs$SocialIdeology)
MPs$EconIdeo01 <- range01(MPs$EconIdeology, na.rm = TRUE)
summary(MPs$EconIdeo01)
MPs$SocioIdeo01 <- range01(MPs$SocialIdeology, na.rm = TRUE)
summary(MPs$SocioIdeo01)
MPs$Ideology01 <- (MPs$SocioIdeo01 + MPs$EconIdeo01)/2
psych::describe(MPs$Ideology01)

#### RECODE: Tenure ####
table(MPs$Tenure1)
MPs$Tenure1 <- factor(MPs$Tenure,
                      levels = c("OneYr",
                                 "FiveYr", "TenYr", "TenYr+"),
                      ordered = TRUE)
MPs$Tenure1 <- reorder(MPs$Tenure1, ref = "OneYr")

#### Please note that sex and age have been redacted to protect participant anonymity as per ethical approval. 
#### Neither variable impacts the inferential analysis in this script.

#### RECODE: PVQ-24 trustworthiness ####
## Other-to-self perceptions ##

# A. Cognitive 

# Ability trust: The public think you are capable of performing your job competently.
table(MPs$Q56_1)
MPs$CogAbT1 <- factor(MPs$Q56_1,
                      levels = c("Strongly Disagree",
                                 "Disagree",
                                 "Slightly Disagree",
                                 "Neither Agree nor Disagree",
                                 "Slightly Agree",
                                 "Agree",
                                 "Strongly Agree"),
                      ordered = TRUE)
levels(MPs$CogAbT1)
table(MPs$CogAbT1)

# Ability distrust: The public think you waste a lot of public money.
table(MPs$Q56_2)
MPs$CogAbD1 <- factor(MPs$Q56_2,
                      levels = c("Strongly Disagree",
                                 "Disagree",
                                 "Slightly Disagree",
                                 "Neither Agree nor Disagree",
                                 "Slightly Agree",
                                 "Agree",
                                 "Strongly Agree"),
                      ordered = TRUE)
levels(MPs$CogAbD1)
table(MPs$CogAbD1)

# Ability trust: The public think you are good at getting the job done.
table(MPs$Q56_3)
MPs$CogAbT2 <- factor(MPs$Q56_3,
                      levels = c("Strongly Disagree",
                                 "Disagree",
                                 "Slightly Disagree",
                                 "Neither Agree nor Disagree",
                                 "Slightly Agree",
                                 "Agree",
                                 "Strongly Agree"),
                      ordered = TRUE)
levels(MPs$CogAbT2)
table(MPs$CogAbT2)

# Ability distrust: The public think you lack technical expertise.
table(MPs$Q56_4)
MPs$CogAbD2 <- factor(MPs$Q56_4,
                      levels = c("Strongly Disagree",
                                 "Disagree",
                                 "Slightly Disagree",
                                 "Neither Agree nor Disagree",
                                 "Slightly Agree",
                                 "Agree",
                                 "Strongly Agree"),
                      ordered = TRUE)
levels(MPs$CogAbD2)
table(MPs$CogAbD2)

# Benevolence trust: The public think you care about people like them.
table(MPs$Q56_5)
MPs$CogBenT1 <- factor(MPs$Q56_5,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogBenT1)
table(MPs$CogBenT1)

# Benevolence distrust: The public think you tend to look after your own interests rather than trying to help others.
table(MPs$Q56_6)
MPs$CogBenD1 <- factor(MPs$Q56_6,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogBenD1)
table(MPs$CogBenD1)

# Benevolence trust: The public think you treat all groups in society fairly.
table(MPs$Q56_7)
MPs$CogBenT2 <- factor(MPs$Q56_7,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogBenT2)
table(MPs$CogBenT2)

# Benevolence distrust: The public think you dont really understand the problems facing ordinary people.
table(MPs$Q56_8)
MPs$CogBenD2 <- factor(MPs$Q56_8,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogBenD2)
table(MPs$CogBenD2)

# Integrity trust: The public think you adhere to a strong moral code.
table(MPs$Q56_9)
MPs$CogIntT1 <- factor(MPs$Q56_9,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogIntT1)
table(MPs$CogIntT1)

# Integrity distrust: The public think you distort the facts to make policies look good.
table(MPs$Q56_10)
MPs$CogIntD1 <- factor(MPs$Q56_10,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogIntD1)
table(MPs$CogIntD1)

# Integrity trust: The public think you tell them the truth.
table(MPs$Q56_11)
MPs$CogIntT2 <- factor(MPs$Q56_11,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogIntT2)
table(MPs$CogIntT2)

# Integrity distrust: The public think you are happy to make promises at elections, but then forget them afterwards.
table(MPs$Q56_12)
MPs$CogIntD2 <- factor(MPs$Q56_12,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$CogIntD2)
table(MPs$CogIntD2)

# B. Affective 

# Positive: The public feel hopeful that you can improve their lives.

table(MPs$Q56_13)
MPs$AffPosT1 <- factor(MPs$Q56_13,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$AffPosT1)
table(MPs$AffPosT1)

# Positive: The public feel confident that you will treat them fairly.
table(MPs$Q56_14)
MPs$AffPosT2 <- factor(MPs$Q56_14,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$AffPosT2)
table(MPs$AffPosT2)

# Positive: The public feel assured that you are doing your best to represent their interests.
table(MPs$Q56_15)
MPs$AffPosT3 <- factor(MPs$Q56_15,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$AffPosT3)
table(MPs$AffPosT3)

# Negative: The public fear that you try to take advantage of them.
table(MPs$Q56_16)
MPs$AffNegD1 <- factor(MPs$Q56_16,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$AffNegD1)
table(MPs$AffNegD1)

# Negative: The public are sceptical of what you say to them.
table(MPs$Q56_17)
MPs$AffNegD2 <- factor(MPs$Q56_17,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$AffNegD2)
table(MPs$AffNegD2)

# Negative: The public feel angry that you show them a lack of respect.
table(MPs$Q56_18)
MPs$AffNegD3 <- factor(MPs$Q56_18,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$AffNegD3)
table(MPs$AffNegD3)

# C. Behavioural intentional 

# Supportive: The public speak openly with you.
table(MPs$Q56_19)
MPs$BehSupT1 <- factor(MPs$Q56_19,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$BehSupT1)
table(MPs$BehSupT1)

# Supportive: The public would vote for you again if they'd had the chance to do so at the last election.
table(MPs$Q56_20)
MPs$BehSupT2 <- factor(MPs$Q56_20,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$BehSupT2)
table(MPs$BehSupT2)

# Supportive: The public seek your help when they are in trouble.
table(MPs$Q56_21)
MPs$BehSupT3 <- factor(MPs$Q56_21,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$BehSupT3)
table(MPs$BehSupT3)

# Monitory: The public monitor your behaviour closely.
table(MPs$Q56_22)
MPs$BehMonD1 <- factor(MPs$Q56_22,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$BehMonD1)
table(MPs$BehMonD1)

# Monitory: The public check whether you have met your electoral promises.
table(MPs$Q56_23)
MPs$BehMonD2 <- factor(MPs$Q56_23,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$BehMonD2)
table(MPs$BehMonD2)

# Monitory: The public double-check what you tell them for misleading information.
table(MPs$Q56_24)
MPs$BehMonD3 <- factor(MPs$Q56_24,
                       levels = c("Strongly Disagree",
                                  "Disagree",
                                  "Slightly Disagree",
                                  "Neither Agree nor Disagree",
                                  "Slightly Agree",
                                  "Agree",
                                  "Strongly Agree"),
                       ordered = TRUE)
levels(MPs$BehMonD3)
table(MPs$BehMonD3)

#### RECODE: Experiments/Treatments ####

#### Experiment 1 - Pandemics (related questions: Q24, Q26) ####
MPs$treat1 <- ifelse(MPs$Q24 != "", "Gain", 
                     ifelse(MPs$Q26 != "", "Loss", NA))
table(MPs$treat1)

## Binary DV (which way do you vote?)
MPs$pandemic1 <- ifelse(MPs$treat1 == "Gain", MPs$Q24, 
                        ifelse(MPs$treat1 == "Loss", MPs$Q26, NA))
table(MPs$pandemic1)
MPs$pandemic1[MPs$pandemic1=="Vote against:program B is implemented"] <- "Vote against (take the risky option)"
MPs$pandemic1[MPs$pandemic1=="Vote against: program B is implemented."] <- "Vote against (take the risky option)"
MPs$pandemic1[MPs$pandemic1=="Vote in favour: program A is implemented"] <- "Vote in favour (avoid the risky option)"
MPs$pandemic1[MPs$pandemic1=="Vote in favour: program A is implemented."] <- "Vote in favour (avoid the risky option)"

MPs$pan01 <- ifelse(MPs$pandemic1 == "Vote against (take the risky option)", 1, 0)
table(MPs$pan01)

#### Experiment 2 - Speaking about a protest movement (related questions: Q31, Q34) ####
MPs$treat2 <- ifelse(MPs$Q31 != "", "Gain", 
                     ifelse(MPs$Q34 != "", "Loss", NA))
table(MPs$treat2)

## Binary DV (which way do you vote?)
MPs$protest1 <- ifelse(MPs$treat2 == "Gain", MPs$Q31, 
                       ifelse(MPs$treat2 == "Loss", MPs$Q34, NA))
table(MPs$protest1)
MPs$pro01 <- ifelse(MPs$protest1 == "Speak in favour", 1, 0)
table(MPs$pro01)

#### Experiment 3 - Voting on new fiscal measures (related questions: Q45, Q49) ####
MPs$treat3 <- ifelse(MPs$Q45 != "", "Gain", 
                     ifelse(MPs$Q49 != "", "Loss", NA))
table(MPs$treat3)

## Binary DV (which way do you vote?)
MPs$economy1 <- ifelse(MPs$treat3 == "Gain", MPs$Q45, 
                       ifelse(MPs$treat3 == "Loss", MPs$Q49, NA))
table(MPs$economy1)
MPs$economy1[MPs$economy1=="Vote against: the fiscal package is not implemented."] <- "Vote against (avoid the risky option)"
MPs$economy1[MPs$economy1=="Vote in favour: the fiscal package is implemented."] <- "Vote in favour (take the risky option)"

MPs$eco01 <- ifelse(MPs$economy1 == "Vote in favour (take the risky option)", 1, 0)
table(MPs$eco01)

#### Descriptives for Experiment DVs ####
# Binary choice by treatment conditions
mytable <- table(MPs$pan01,MPs$treat1)
prop.table(mytable, 2)*100

mytable <- table(MPs$pro01,MPs$treat2)
prop.table(mytable, 2)*100

mytable <- table(MPs$eco01,MPs$treat3)
prop.table(mytable, 2)*100


#### Fisher's test for treatment effects ####
fisher.test(table(MPs$pan01,MPs$treat1))
fisher.test(table(MPs$pro01,MPs$treat2))
fisher.test(table(MPs$eco01,MPs$treat3))

#### Simulated probabilities for treatment effects ####
library(finalfit)
library(Zelig)

## Experiment one
Logit1 <- glm(pan01 ~ treat1, data = MPs, family = "binomial")
summary(Logit1, cluster = c("Country"))
newdata <- with(MPs, data.frame(treat1 = c("Gain","Loss")))
boot_predict(Logit1, newdata, R=1000, boot_compare = FALSE,
             digits = c(2,3))
## 1000 simulated iterations suggests that 58% of politicians in the loss treatment take the risk
## compared to just 12% in the gain treatment.

## Use Zelig to plot range of predicted probabilities in each condition.
# re-estimate model
z5_1 <- zelig(pan01 ~ as.factor(treat1), model = "logit", data = MPs, cite = FALSE)
# model summary
summary(z5_1)
# set treat1 to 1 and 0
z5_1 <- setx(z5_1, treat1 = "Gain")
z5_1 <- setx1(z5_1, treat1 = "Loss")
# run simulations and estimate quantities of interest
z5_1 <- sim(z5_1, num = 1000)
# model summary
summary(z5_1)
# plot values summary
plot(z5_1)

# extract predicted probabilities
Pred <- z5_1$get_qi(qi ="ev", xvalue = "x")
Pred1 <- z5_1$get_qi(qi ="ev", xvalue = "x1")
Pred <- as.data.frame(Pred)
Pred$Treat <- "Gain"
Pred1 <- as.data.frame(Pred1)
Pred1$Treat <- "Loss"
Preds <- rbind(Pred, Pred1) ## use this for visualisation (below)

## Experiment two
Logit1 <- glm(pro01 ~ treat2, data = MPs, family = "binomial")
summary(Logit1, cluster = c("Country"))
newdata <- with(MPs, data.frame(treat2 = c("Gain","Loss")))
boot_predict(Logit1, newdata, R=1000, boot_compare = FALSE,
             digits = c(2,3))
## 1000 simulated iterations suggests that 55% of politicians in the loss treatment take the risk
## compared to 47% in the gain treatment.

## Use Zelig to plot range of predicted probabilities in each condition.
# re-estimate model
z5_1 <- zelig(pro01 ~ as.factor(treat2), model = "logit", data = MPs, cite = FALSE)
# model summary
summary(z5_1)
# set treat1 to 1 and 0
z5_1 <- setx(z5_1, treat2 = "Gain")
z5_1 <- setx1(z5_1, treat2 = "Loss")
# run simulations and estimate quantities of interest
z5_1 <- Zelig::sim(z5_1, num = 1000)
# model summary
summary(z5_1)
# plot values summary
plot(z5_1)

# extract predicted probabilities
Pred <- z5_1$get_qi(qi ="ev", xvalue = "x")
Pred1 <- z5_1$get_qi(qi ="ev", xvalue = "x1")
Pred <- as.data.frame(Pred)
Pred$Treat <- "Gain"
Pred1 <- as.data.frame(Pred1)
Pred1$Treat <- "Loss"
Preds1 <- rbind(Pred, Pred1) ## use this for visualisation (below)

## Experiment three
Logit1 <- glm(eco01 ~ treat3, data = MPs, family = "binomial")
summary(Logit1, cluster = c("Country"))
newdata <- with(MPs, data.frame(treat3 = c("Gain","Loss")))
boot_predict(Logit1, newdata, R=1000, boot_compare = FALSE,
             digits = c(2,3))
## 1000 simulated iterations suggests that 38% of politicians in the loss treatment take the risk
## compared to 33% in the gain treatment.

## Use Zelig to plot range of predicted probabilities in each condition.
# re-estimate model
z5_1 <- zelig(eco01 ~ as.factor(treat3), model = "logit", data = MPs, cite = FALSE)
# model summary
summary(z5_1)
# set treat1 to 1 and 0
z5_1 <- setx(z5_1, treat3 = "Gain")
z5_1 <- setx1(z5_1, treat3 = "Loss")
# run simulations and estimate quantities of interest
z5_1 <- Zelig::sim(z5_1, num = 1000)
# model summary
summary(z5_1)
# plot values summary
plot(z5_1)

# extract predicted probabilities
Pred <- z5_1$get_qi(qi ="ev", xvalue = "x")
Pred1 <- z5_1$get_qi(qi ="ev", xvalue = "x1")
Pred <- as.data.frame(Pred)
Pred$Treat <- "Gain"
Pred1 <- as.data.frame(Pred1)
Pred1$Treat <- "Loss"
Preds2 <- rbind(Pred, Pred1) ## use this for visualisation (below)

## Visual summaries of treatment effects

Preds$Exp <- "'Asian' Disease"
Preds1$Exp <- "PR Conundrum"
Preds2$Exp <- "Economic Turmoil"
PredProb <- rbind(Preds, Preds1, Preds2)

PredProb %>%
  ggplot(aes(x=V1, y=..scaled.., color=Treat, linetype = Treat)) +
  geom_density(alpha = 0.6, size = 2.5, colour = "black") +
  theme_bw() +
  facet_wrap(~Exp) +
  labs(x = "Probability of risk taking (%)",
       y = "Density",
       caption = "Density plots of simulated predicted probabilities for each treatment condition. \nCalculated from univariate logistic regression with treatment conditions as the only independent variable.\n(Iterations = 1,000 per panel.)") +
  theme(legend.position = "right",
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        axis.title = element_text(size = 18, face = "bold"),
        axis.text.x = element_text(size = 18),
        axis.text.y = element_text(size = 18),
        plot.caption = element_text(size = 12),
        strip.text = element_text(size = 12, face = "italic")) 

#### Treatment check: qualitative experiment reflections ####

MPs$Q69[MPs$Q69 == ""] <- NA
table(is.na(MPs$Q69)) ## 91 participants provided qualitative reflections on the experiments.

## install appropriate packages
library(readtext)
library(quanteda)
library(quanteda.textstats)
library(tm)

# Note: export qualitative answers to text file (Q69 in the dataset)
filepath <- "INSERT APPROPRIATE FILE DIRECTORY"
rt <- readtext(filepath, text_field = "texts")
rt

full_corpus <- corpus(rt) # create quanteda corpus
summary(full_corpus) # 70 sentences comprising 1719 tokens

dtm <- dfm(full_corpus, tolower = TRUE, remove_numbers = TRUE, 
           remove_punct = TRUE, remove_separators = TRUE,
           stem = FALSE,remove = stopwords("english"))
dtm

# plot 20 most frequent words (Overall frequencies)
library("ggplot2")
freq <- textstat_frequency(dtm)
head(freq, 10)
ggplot(freq[1:20, ], aes(x = reorder(feature, frequency), y = frequency)) +
  geom_point(shape = 18, colour = "black", fill = "black", size = 3) + 
  ylab ('Total Frequency') +
  xlab ('Frequent Terms') +
  ggtitle('Thinking about your answers so far, what influenced your decisions?') +
  theme(plot.title = element_text(size=16, face='bold'),
        axis.text=element_text(size=16),
        axis.title=element_text(size=16),
        panel.background = element_rect(fill = 'white', colour = 'black')) +
  coord_flip() 

## Plot relative frequencies of terms across all participant answers
dfm_weight <- full_corpus %>% 
  dfm(tolower = TRUE, remove_numbers = TRUE, 
      remove_punct = TRUE, remove_separators = TRUE,
      stem = FALSE, remove = stopwords("english"), remove_url = TRUE, remove_twitter = TRUE) %>% 
  dfm_group(groups = "Document") %>%
  dfm_weight(scheme = "prop")
freq_weight <- textstat_frequency(dfm_weight, n = 1,
                                  groups = "Document")

## discounting suggestable stopwords, the most common terms are:
# 1. parties/party
# 2. decision
# 3. whip
# 4. outcome
# 5. risk
# 6. people

## search for these words in context
## Locate phrases using response (revise and update accordingly for multiple searches)
head(kwic(full_corpus, c("part*"), window = 10, valuetype = "glob"))

## Look back in the datafile for match examples with participants
# Three examples from each country:

# Canadian MP: The first one, no easy answer. The second one, I couldn't return to the starting point of the question, so assumed that the protest movement was positive (so party discipline matters less to me), and on third question since all else was equal, the party's view tipped balance.
# Canadian MP: 1. We chose to carry the burden of decision making. 2. We can disagree and debate on an issue in the caucus, but outside of caucus room we are united. 3. Polling plays a smaller part in my decision making. 
# Canadian MP: For the sickness, I don't like to gamble with peoples lives. On issues, it's better to take a stand if you can than dodge it. For fiscal measures, I don't like stimulus or bailout packages generally, so I'd just as soon follow party lines.
# UK MP: Q1: Probability of greater deaths; Q2: need to make an informed decision dependent on the subject; Q3: it is a whipped vote!
# UK MP: Q1: Chance to save all people needs to be enlarged; Q2: you have to lead; Q3: the vote is whipped
# UK MP: Q1: I voted for the option which holds out hope that nobody may die, even though the likely number of deaths would be higher; Q2: the expected outcome of each alternative is no change to party policy. I decided it was worth risking the choice with a potential policy upside; Q3: the expected outcome of the package on the economy is neutral, so I am happy to comply with the whip and vote against the package.
# SA MP: On free vote I will try and follow my own value system. On whipped decisions, there is no choice, party discipline must be maintained. 
# SA MP: 1. more chance to save lives and bought time to increase odds 2. Spoke for itself 3. I represent party and vote with it. 
# SA MP: The first question: if we can save 2000 people now, whose to say there will not be a break through later where we can save the other 3000. Perhaps if there was a 50% chance opposed to 35% chance to save everyone, I would opt for that. the Second question: I need to trust that the decision to not support was discussed throughly in our caucus with all questions posed, and therefore support the decision of the caucus unless it posed an ethical risk.

#### Processing/Robustness checks on PVQ-24 data ####

#### Ordinary scale scores for Trust, distrust, mistrust and alphas ####
## Create generic trust, distrust and mistrust variables
MPs$TRUST <- (as.numeric(MPs$CogAbT1) + as.numeric(MPs$CogAbT2) +
                as.numeric(MPs$CogBenT1) + as.numeric(MPs$CogBenT2) +
                as.numeric(MPs$CogIntT1) + as.numeric(MPs$CogIntT2)  +
                as.numeric(MPs$AffPosT1) + as.numeric(MPs$AffPosT2) +
                as.numeric(MPs$AffPosT3) + as.numeric(MPs$BehSupT1) + as.numeric(MPs$BehSupT2) +
                as.numeric(MPs$BehSupT3))/12
psych::describe(MPs$TRUST)
Trust <- data.frame(as.numeric(MPs$CogAbT1) , as.numeric(MPs$CogAbT2) ,
                    as.numeric(MPs$CogBenT1) , as.numeric(MPs$CogBenT2) ,
                    as.numeric(MPs$CogIntT1) , as.numeric(MPs$CogIntT2)  ,
                    as.numeric(MPs$AffPosT1) , as.numeric(MPs$AffPosT2) ,
                    as.numeric(MPs$AffPosT3) , as.numeric(MPs$BehSupT1) , as.numeric(MPs$BehSupT2) ,
                    as.numeric(MPs$BehSupT3))
alpha(Trust, na.rm = TRUE) # .87

MPs$DISTRUST <- (as.numeric(MPs$CogAbD1) + as.numeric(MPs$CogAbD2) +
                   as.numeric(MPs$CogBenD1) + as.numeric(MPs$CogBenD2) +
                   as.numeric(MPs$CogIntD1) + as.numeric(MPs$CogIntD2) +
                   as.numeric(MPs$AffNegD1) + as.numeric(MPs$AffNegD2) +
                   as.numeric(MPs$AffNegD3))/9
Distrust <- data.frame(as.numeric(MPs$CogAbD1) , as.numeric(MPs$CogAbD2) ,
                       as.numeric(MPs$CogBenD1) , as.numeric(MPs$CogBenD2) ,
                       as.numeric(MPs$CogIntD1) , as.numeric(MPs$CogIntD2) ,
                       as.numeric(MPs$AffNegD1) , as.numeric(MPs$AffNegD2) ,
                       as.numeric(MPs$AffNegD3))
alpha(Distrust, na.rm = TRUE) # .9

MPs$MISTRUST <- (as.numeric(MPs$BehMonD1) + as.numeric(MPs$BehMonD2) +
                   as.numeric(MPs$BehMonD3))/3
Mistrust <- data.frame(as.numeric(MPs$BehMonD1) , as.numeric(MPs$BehMonD2) ,
                       as.numeric(MPs$BehMonD3))
alpha(Mistrust) # .55 Lower alpha but expected with only three items


#### Step one: MDS analysis ####
# a. isolate desired variables in new dataset with only positive values for each variable:
MDSPerception <- data.frame(MPs$CogAbT1, MPs$CogAbT2,
                            MPs$CogBenT1, MPs$CogBenT2,
                            MPs$CogIntT1, MPs$CogIntT2,
                            MPs$CogAbD1, MPs$CogAbD2,
                            MPs$CogBenD1, MPs$CogBenD2,
                            MPs$CogIntD1, MPs$CogIntD2,
                            MPs$AffPosT1, MPs$AffPosT2,
                            MPs$AffPosT3, MPs$AffNegD1,
                            MPs$AffNegD2, MPs$AffNegD3,
                            MPs$BehSupT1, MPs$BehSupT2,
                            MPs$BehSupT3, MPs$BehMonD1,
                            MPs$BehMonD2, MPs$BehMonD3)

colnames(MDSPerception) <- c("CAT1", "CAT2", "CBT1", "CBT2", "CIT1", "CIT2",
                             "CAD1", "CAD2", "CBD1", "CBD2", "CID1", "CID2",
                             "APT1", "APT2", "APT3", "AND1", "AND2", "AND3",
                             "BST1", "BST2", "BST3", "BMD1", "BMD2", "BMD3")

# b. Load vegan package and calculate MDS package:
install.packages("vegan")
library("vegan")
data("MDSPerception")

# the following code:
# a. removes transformations for ecological data built into 'vegan'
# b. specifies the number of dimensions (k)
# c. requests number of attempts to reduce 'stress' in the ordination space (limit this to 30)
# ... stress refers to the variablility in the dissimilarity matrix (using all data) explained by the ordination distances
# ... a value of up to .2 is generally acceptable and indicates 80% variability explained

# I want to examine 3 dimensions - i.e. theoretical expectations of trust, distrust, and mistrust
MDSPerception <- na.omit(MDSPerception)

Perception.data <- metaMDS(MDSPerception, distance='euclidean', k=3, 
                           trymax=30, autotransform=FALSE, 
                           wasscores=TRUE, noshare=FALSE)

# c. Visualise the stress between ordination and dissimilarity matrices on a stressplot
# Two stress values are provided:
# Non-metric fit - this is the modified Kruskal's stress value and is 1−R2‾‾‾‾‾‾‾√
# Linear fit - this is the R2 value of the correlation between the ordination values and 
# the ordination values predicted from the (monotonic) regression line.
# ... looking for 
stressplot(Perception.data)

# d. plot the ordination space with ggplot and customise!
install.packages("remotes")
remotes::install_github("jfq3/ggordiplots")
library(ggordiplots)

# Start by running base MDS plot
pl <- ordiplot(Perception.data, type = "none", xlab = "Dimension 1", ylab = "Dimension 2")
# Then extract species data and turn it to a dataframe
# Add an identifying variable (i.e. item names and group names)
# Plot with ggplot and add shapes, colours, labels (i.e. variability explained by each dimension)

#1. 
ord.plot <- pl$species
ord.plot <- as.data.frame(ord.plot)
# 2. 
item.names <- c("CAT1", "CAT2", "CBT1", "CBT2", "CIT1", "CIT2",
                "CAD1", "CAD2", "CBD1", "CBD2", "CID1", "CID2",
                "APT1", "APT2", "APT3", "AND1", "AND2", "AND3",
                "BST1", "BST2", "BST3", "BMD1", "BMD2", "BMD3")
type.trust <- c("Trust", "Trust", "Trust", "Trust", "Trust", "Trust",
                "Distrust", "Distrust", "Distrust", "Distrust","Distrust", "Distrust",
                "Trust", "Trust", "Trust", "Distrust", "Distrust", "Distrust",
                "Trust", "Trust", "Trust", "Distrust", "Distrust", "Distrust")
type.facet <- c("Cognitive", "Cognitive", "Cognitive", "Cognitive", "Cognitive", "Cognitive",
                "Cognitive", "Cognitive", "Cognitive", "Cognitive", "Cognitive", "Cognitive",
                "Affective", "Affective", "Affective", "Affective", "Affective", "Affective",
                "Behavioural-intentional", "Behavioural-intentional", "Behavioural-intentional",
                "Behavioural-intentional", "Behavioural-intentional", "Behavioural-intentional")
# 3.
ord.plot$items <- item.names
ord.plot$trust <- type.trust
ord.plot$facet <- type.facet
#4.
# Plot based on three dimensional MDS:
Perception.data$species

p <- ggplot(data = ord.plot) +
  aes(x = NMDS1, y = NMDS2, shape = trust, colour = facet) + 
  geom_point(size = 8) + 
  labs(x = "Dimension (axis) 1", y = "Dimension (axis) 2", 
       #title = "Multi-dimensional scaling of PTB-24 responses",
       #subtitle = "Politicians' perceptions of public trust judgements",
       caption = "monoMDS \nDistance = euclidean \nDimensions = 3, \nStress(30 computations) = 0.116") +
  scale_color_manual(values = c("#BDBDBD", "#969696", "#252525")) +
  guides(color = guide_legend("Judgement Type")) +
  guides(shape = guide_legend("Intended Latent Concept")) +
  theme_bw() +
  theme(
    #axis.title.y = element_text(vjust= 2.0),
    #axis.title.x = element_text(vjust= -1.0),
    axis.title = element_text(size = 16, face = "bold"),
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 12, face = "bold"), 
    plot.caption = element_text(size = 12),
    plot.title = element_text(size = 16, face = "bold", vjust = 2.0),
    plot.subtitle = element_text(size = 12, face = "bold", vjust = 1.0)) 
p

library(RColorBrewer)
brewer.pal(n=8, name = "Greys")

## Add text labels?
p + geom_text(aes(label=rownames(ord.plot)),hjust=0, vjust=0)

install.packages("ggrepel")
library(ggrepel)
p + geom_label_repel(aes(label = rownames(ord.plot)),
                     box.padding   = 0.35, 
                     point.padding = 0.5,
                     segment.color = 'grey50') 

#### Step two: PCA analysis ####
## For MPs only 
which( colnames(MPs)=="CogAbT1" ) # check column numbers of first and last items in PTB
which( colnames(MPs)=="BehMonD3" )

Fit <- principal(MPs[151:174], nfactors = 3, rotate = "promax", missing = TRUE, impute = "median")
print(Fit, digits=2, cutoff=.3) # Good fit but some weak loadings
# The following loaded below .2/.3 - remove and re-run?
which( colnames(MPs)=="CogBenT1" )
which( colnames(MPs)=="CogBenT2" )
which( colnames(MPs)=="CogIntT1" )
which( colnames(MPs)=="CogIntT2" )
which( colnames(MPs)=="AffPosT2" )
MPs1 <- MPs[ -c(155,157,159,161,164) ]

which( colnames(MPs1)=="CogAbT1" )
which( colnames(MPs1)=="BehMonD3" )

Fit <- principal(MPs1[151:169], nfactors = 3, rotate = "promax", missing = TRUE, impute = "median")
print(Fit, digits=2, cutoff=.3)
Fact <- data.frame(Fit$scores)

# Add factor scores to dataframe
MPs$TRUSTFACT <- Fact$RC2
MPs$DISTRUSTFACT <- Fact$RC1
MPs$MISTRUSTFACT <- Fact$RC3
# Rescale 0-1
MPs$TRUSTFACT01 <- range01(MPs$TRUSTFACT, na.rm = TRUE)
MPs$DISTRUSTFACT01 <- range01(MPs$DISTRUSTFACT, na.rm = TRUE)
MPs$MISTRUSTFACT01 <- range01(MPs$MISTRUSTFACT, na.rm = TRUE)

# Check descriptives and distributions
psych::describe(MPs$TRUSTFACT)
ggplot(MPs) +
  aes(x = TRUSTFACT) +
  geom_histogram()

psych::describe(MPs$DISTRUSTFACT)
ggplot(MPs) +
  aes(x = DISTRUSTFACT) +
  geom_histogram()

psych::describe(MPs$MISTRUSTFACT)
ggplot(MPs) +
  aes(x = MISTRUSTFACT) +
  geom_histogram()

#### Step three: CFA with single trust, distrust and mistrust dimensions ####
MODEL <- 
  'TRUSTCFA =~ CogAbT1 + CogAbT2 + CogBenT1 + CogBenT2 + CogIntT1 + CogIntT2 +
AffPosT1 + AffPosT2 + AffPosT3 + BehSupT1 + BehSupT2 + BehSupT3
DISTRUSTCFA =~ CogAbD1 + CogAbD2 + CogBenD1 + CogBenD2 + CogIntD1 + CogIntD2 +
AffNegD1 + AffNegD2 + AffNegD3 
MISTRUSTCFA =~ BehMonD1 + BehMonD2 + BehMonD3'

Model3 <- cfa(MODEL, data = MPs, missing = "fiml", estimator = "MLR", cluster = "Country")
summary(Model3, fit.measures = TRUE, standardized = TRUE)

## Remove weak predictors - items loading less than .3 in PCA
MODEL <- 
  'TRUSTCFA =~ CogAbT1 + CogAbT2 + CogBenT1 + CogBenT2 + CogIntT1 + CogIntT2 +
AffPosT1 + AffPosT2 + AffPosT3 + BehSupT2 + BehSupT3
DISTRUSTCFA =~ CogAbD2 + CogBenD2 + CogIntD1 + CogIntD2 +
AffNegD1 + AffNegD2  
MISTRUSTCFA =~ BehMonD1 + BehMonD2 + BehMonD3'

Model3 <- cfa(MODEL, data = MPs, missing = "fiml", estimator = "MLR", cluster = "Country")
summary(Model3, fit.measures = TRUE, standardized = TRUE)

#### Mixed Effects Regression ####
library(sjstats) # for ICCs
library(MuMIn) # for R2
install.packages("optimx")
library(optimx)

#### 1. Experiment One: Health Crisis ####
# Country random intercepts only
Logit1 <- glmer(pan01 ~ (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit1)
r.squaredGLMM(Logit1) 

# + trust variables and treatment effects
Logit12 <- glmer(pan01 ~ treat1 + TRUSTFACT + DISTRUSTFACT + MISTRUSTFACT + (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit12)
r.squaredGLMM(Logit12) 

# + interactions
Logit13 <- glmer(pan01 ~ treat1*TRUSTFACT + treat1*DISTRUSTFACT + treat1*MISTRUSTFACT + (1 |Country), 
                 data = MPs, family = "binomial", weights = weight, glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(Logit13)
r.squaredGLMM(Logit13) 

# + random slopes and controls
Logit14 <- glmer(pan01 ~ treat1*TRUSTFACT + treat1*DISTRUSTFACT + treat1*MISTRUSTFACT +
                   Safeseat + Ideology01 + Tenure1 + 
                   (1 + TRUSTFACT + DISTRUSTFACT + MISTRUSTFACT |Country), 
                 data = MPs, family = "binomial", weights = weight, 
                 glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(Logit14)
r.squaredGLMM(Logit14) 
cc <- coef(summary(Logit14))
cc <- cbind(cc,adjust.p=p.adjust(cc[,"Pr(>|z|)"],"bonferroni")) ## bonferroni
## note that even after corrections, trust:treatment interaction remains sig at p<.10

#### 2. Experiment Two: PR Conundrum ####
# Country random intercepts only
Logit2 <- glmer(pro01 ~ (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit2)
r.squaredGLMM(Logit2) 

# + trust variables and treatment effects
Logit21 <- glmer(pro01 ~ treat2 + TRUSTFACT + DISTRUSTFACT + MISTRUSTFACT + (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit21) # no meaningful main effects
r.squaredGLMM(Logit21) 

# + interactions
Logit22 <- glmer(pro01 ~ treat2*TRUSTFACT + treat2*DISTRUSTFACT + treat2*MISTRUSTFACT + (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit22) # only treatments significant at p<.10
r.squaredGLMM(Logit22) 

# + random slopes and controls
Logit23 <- glmer(pro01 ~ treat2*TRUSTFACT + treat2*DISTRUSTFACT + treat2*MISTRUSTFACT +
                   Safeseat + Ideology01 + Tenure1 + 
                   (1 + TRUSTFACT + DISTRUSTFACT + MISTRUSTFACT|Country), 
                 data = MPs, family = "binomial", weights = weight, 
                 glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(Logit23)
r.squaredGLMM(Logit23)

#### 3. Experiment Three: PR Conundrum ####
# Country random intercepts only
Logit3 <- glmer(eco01 ~ (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit3)
r.squaredGLMM(Logit3) 

# + trust variables and treatment effects
Logit31 <- glmer(eco01 ~ treat3 + TRUSTFACT + DISTRUSTFACT + MISTRUSTFACT + (1 |Country), data = MPs, family = "binomial", weights = weight)
summary(Logit31) # clear treatment effect
r.squaredGLMM(Logit31) 

# + interactions
Logit32 <- glmer(eco01 ~ treat3*TRUSTFACT + treat3*DISTRUSTFACT + treat3*MISTRUSTFACT + (1 |Country), 
                 data = MPs, family = "binomial", weights = weight, 
                 glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(Logit32) # strong distrust:treatment interaction
r.squaredGLMM(Logit32) 
cc <- coef(summary(Logit32))
cc <- cbind(cc,adjust.p=p.adjust(cc[,"Pr(>|z|)"],"bonferroni")) ## bonferroni
## still significant at p<.05 after correction

# + random slopes and controls
Logit33 <- glmer(eco01 ~ treat3*TRUSTFACT + treat3*DISTRUSTFACT + treat3*MISTRUSTFACT +
                   Safeseat + Ideology01 + Tenure1 + 
                   (1 + TRUSTFACT + DISTRUSTFACT + MISTRUSTFACT|Country), 
                 data = MPs, family = "binomial", weights = weight, 
                 glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
summary(Logit33)
r.squaredGLMM(Logit33) 
cc <- coef(summary(Logit33))
cc <- cbind(cc,adjust.p=p.adjust(cc[,"Pr(>|z|)"],"bonferroni")) ## bonferroni
## distrust:treatment interaction drops to p<.15 in complex model
## BUT small sample size may account for this given Bonferroni can be overly harsh

#### Regression tables for paper ####
library(sjPlot)

# 1. Experiment One
tab_model(Logit1, Logit12, Logit13, Logit14, show.est = TRUE)

# 2. Experiment Two
tab_model(Logit2, Logit21, Logit22, Logit23, show.est = TRUE)

# 3. Experiment Three
tab_model(Logit3, Logit31, Logit32, Logit33, show.est = TRUE)

# 4. Complex models together
tab_model(Logit14, Logit23, Logit33, show.est = TRUE)


